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Numerical treatment of the Filament Based Lamellipodium 

Model (FBLM) 

A. Manhart * * * , D. Oelz \ C. Schmeiser +, N. Sfakianakis § , 


Abstract 

We describe in this work the numerical treatment of the Filament Based Lamellipodium Model (FBLM). 
The model itself is a two-phase two-dimensional continuum model, describing the dynamics of two interacting 
families of locally parallel F-actin filaments. It includes, among others, the bending stiffness of the filaments, 
adhesion to the substrate, and the cross-links connecting the two families. The numerical method proposed 
is a Finite Element Method (FEM) developed specifically for the needs of these problem. It is comprised of 
composite Lagrange-Hermite two dimensional elements defined over two dimensional space. We present some 
elements of the FEM and emphasise in the numerical treatment of the more complex terms. We also present 
novel numerical simulations and compare to in-vitro experiments of moving cells. 


1 Introduction 


The lamellipodium is a flat cell protrusion mechanism functioning as a motility organelle in protrusive cell mi¬ 
gration [31]. It is a very dynamic stmcture mainly consisting of a network of branched actin filaments. These are 
semi-elastic rods that represent the polymer form of the protein actin. They are continuously remodeled by poly¬ 
merization and depolymerization and therefore undergo treadmilling [2], Actin associated cross-linker proteins 
and myosin motor proteins integrate them into the lamellipodial meshwork which plays a key role in cell shape 
stabilization and in cell migration. Different modes of cell migration result from the interplay of protrusive forces 
due to polymerization, actomyosin dependent contractile forces and regulation of cell adhesion [11]. 

The first modeling attempts have resolved the interplay of protrusion at the front and retraction at the rear in a 
one-dimensional spatial setting [1,7], Two-dimensional continuum models were developed in order to include the 
lateral flow of F-actin along the leading edge of the cell into the quantitative picture. Those models can explain 
characteristic shapes of amoeboid cell migration [25, 24] on two-dimensional surfaces as well as the transition to 
mesenchymal migration [26], 

One of the still unresolved scientific questions concerns the interplay between macroscopic observables of cell 
migration and the microstructure of the lamellipodium meshwork. Specialised models had been developed sepa¬ 
rately from the continuum models to track microscopic information on filament directions and branching structure 
[13, 27, 10]. However, solving fluid-type models that describe the whole cytoplasm while retaining some informa¬ 
tion on the microstructure of the meshwork has turned out to be challenging. One approach is to formulate hybrid 
models [17], another one to directly formulate models on the computational, discrete level [16], Recently the 
approach to directly formulate a computational model has been even extended into the three-dimensional setting 
making use of a finite element discretization [18]. 

In an attempt to create a simulation framework that addresses the interplay of macroscopic features of cell migra¬ 
tion and the meshwork structure the Filament Based Lamellipodium Model (FBLM) has been developed. It is a 
two-dimensional, two-phase, anisotropic continuum model for the dynamics of the lamellipodium network which 
retains key directional information on the filamenteous substructure of this meshwork [23], 
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Figure 1: Graphical representation of (3); showing here the lamellipodium Q(t) “produced” by the mappings F* and the 
crossing-filament domain C. 


The model has been derived from a microscopic description based on the dynamics and interaction of individual 
filaments [21], and it has by recent extensions [15] reached a certain state of maturity. Since the model can be 
written in the form of a generalized gradient flow, numerical methods based on optimization techniques have 
been developed [22, 23], Numerical efficiency had been a shortcoming of this approach. This has led to the 
development of a Finite Element numerical method which is presented in this article alongside simulations of a 
series of migration assays. 


2 Mathematical modeling 

In this section the FBLM will be sketched. For more detail see [15], The main unknowns of the model are 
the positions of the actin filaments in two locally parallel families (denoted by the superscripts + and -). Each 
of these families covers a topological ring with all individual filaments connecting the inner boundary with the 
outer boundary. The outer boundaries are the physical leading edge and therefore identical, whereas the inner 
boundaries of the two families are artificial and may be different. Filaments are labelled by a e [0, 2n), where the 
interval represents a one-dimensional torus. The maximal arclength of the filaments in an infinitesimal element 
da of the +-family at time t is denoted by L ± (o', t), and an arclength parametrization of the filaments is denoted by 
[F^cr, s, t) : -L ± (a, t) < s < 0} c R 2 , where the leading edge corresponds to s — 0, i.e. 

[FV, 0,0 : 0 < a < 2n] = [F“(or, 0, t) : 0 < a < 2n] ft, (1) 

which together with 

|a j F ± (a,s,f)| = l V(<W), (2) 

constitutes constraints for the unknowns F*. The second constraint is connected to an inextensibility assumption 
on the filaments, which implies that s can also be interpreted as a monomer counter along filaments. 

We expect that different filaments of the same family do not intersect each other, and each plus-filament crosses 
each minus-filament at most once. The first condition is guaranteed by det(3 a .F ± , djF*) > 0, where the sign 
indicates that the labelling with increasing a is in the clockwise direction. The second condition uniquely defines 
s ± = s ± (a + ,a~,t ) such that F + (or + ,s + ,t) = F~(a - , s~, t), for all ( a + ,a~) e C(t), the set of all pairs of crossing 
filaments. As a consequence, there are coordinate transformations t/A : {a + , .s’ 1 ) (a*, s ± ) such that 

F T = F ± o . 
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In the following, we shall concentrate on one of the two families and skip the superscripts except that the other 
family is indicated by the superscript *. The heart of the FBLM is the force balance 

0 = (i^F) - d.s MincxidyF) + /A/D^F + d s (p(p)d a F^ - d a (p(p)d s F^) 

bending in-extensibility adhesion pressure 

±ds (wW - ^P 1 ) + rp 7 y (AF - p;f*) , (3) 

twisting stretching 


where the notation = (7-j, F 2 ) ± = (~F 2 , F\) has been. For fixed ,v and t, the function rj(a, s,t), is the number 
density of filaments of length at least —s at time t with respect to a. Its dynamics and that of the maximal length 
L(a, t) will not be discussed here. It can be modeled by incorporating the effects of polymerization, depolymer¬ 
ization, branching, and capping (see [15]). We only note that, faster polymerization (even locally) leads to wider 
lamellipodia. 

The first term on the right hand side of (3) describes the filaments' resistance against bending with the stiffness pa¬ 
rameter p B > 0. The second term is a tangential tension force, which arises from incorporating the inextensibility 
constraint (2) with the Lagrange multiplier Ti next (a, s, t). The third term describes friction of the filament network 
with the nonmoving substrate (see [21] for its derivation as a macroscopic limit of the dynamics of transient elastic 
adhesion linkages). Since filaments polymerize at the leading edge with the polymerization speed v(a, 1) > 0, they 
are continuously pushed into the cell with that speed, and the material derivative 


Z),F := <3,F - v3 s F 


is the velocity of the actin material relative to the substrate. For the modeling of v see [15]. 

The second line of (3) models a pressure effect caused by Coulomb repulsion between neighboring filaments of 
the same family with pressure pip), where the actin density in physical space is given by 


n 

det(3 ff F, djF) 


(4) 


Finally, the third line of (3) models the interaction between the two families caused by transient elastic cross-links 
and/or branch junctions. The first term describes elastic resistance against changing the angle </> = arccos(5 s F ■ 
diF*) between filaments away from the equilibrium angle </>o of the equilbrium conformation of the cross-linking 
molecule. The last term describes friction between the two families analogously to friction with the substrate. The 
friction coefficients have the form 


pFS = f,s 


da * 
ds 


with p T - s > 0 and the partial derivative refers to the coordinate transformation 1 //, which is also used when 
evaluating partial derivatives of F*. 


The system (3) is considered subject to the boundary conditions 


-p B d s (^“F) - pip)d a F- 1 + 77,li„ext3 s F + pTl*p T {(f> - (/>Q)d s F L 

ft7 (/tan(<*)d s F + / inn (a)V(a)), for s = -L , 
\ ±4tetherV, for S = 0 , 

77<9 x F = 0, for s = -L, 0 . 


(5) 


The terms in the second line are forces applied to the filament ends. The force in the direction v orthogonal to the 
leading edge at s = 0 arises from the constraint (1) with the Lagrange parameter Tether- Its biological interpretation 
is due to tethering of the filament ends to the leading edge. The forces at the inner boundary s = -L are models of 
the contraction effect of actin-myosin interaction in the interior region (see [15] for details). 
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Figure 2: Discretized lamellipodium (left) and lamellipodium fragment (right). 


3 Numerical method 

We present in this section some elements of the Finite Element (FE) method that we have developed for (3). We 
prefer in this work a more computational approach to highlight the implementation difficulties encountered, and to 
emphasize on the “special” treatment that some of the terms necessitate [3, 9]. We refer to [14] for further details 
and a thorough theoretical presentation of the FEM. 

3.1 Discrete formulation 

As previously, we skip the superscripts (+) except for those of the other family that we indicate by *. 

The domain B of F (see also Fig. 1) is discretized into isodynamous computational cells of the form 


Cij = [ai,a i+ 1 ) x [sj, s j+ i) , 


( 6 ) 


where a, v \ - a ; = Acr and ,v,+ \ - Si = As for i = 1... N a — 1, and j = 1... N s - 1. We denote by Tao-.Aj the 
discretization of B. Let now, C represent the generic cell of 7~A a ,As> i.e. 


C = [a\,af) x [si, S 2 ), with a 2 - u\ = A a > 0, S 2 - = As > 0 , 


(7) 


with its vertices in lexicographic order: 


Vi = (ai,si), V 2 = (ai,s 2 ), V 3 = (a 2 , si), V 4 = (a 2 , s 2 ) . 


( 8 ) 


The composite Lcigrcmge-Hennite scalar polynomial space over the cell C is defined as: 



(9) 


where P][[-] denotes the vector space of scalar real polynomials, with real coefficients, of a single variable in the 
corresponding component of C, and degree at most k. 

Remark 3.1. The higher smoothness of the elements of r Vc along the ,v-direction is necessitated by the higher 
order ^-derivatives of F in (3). 


We also set, for ( a , s) e C, the shape functions: 



(10) 


G^(s) = -G^(s 1 + s 2 - s). 
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Figure 3: Graphical representation of the Lagrange-Herrnite shape functions ( 12 ). Each one of the shape functions attains 
the value 1 in one degree of freedom, and 0 on all the rest. 


which satisfy: 


Lf(a 0 = 1, 

TjOi) = 0, 

Gf(o) = 1, 

Lf(a 2 ) = 0, 
Ln(a 2 ) = 1, 
G|(.v 2 ) = 0, 

|G 1 (s 1 ) = 0, 

fG!te)=0, 

G$(s 1 ) = 0, 

G c 2 (s 2 ) = 0, 

|G 2 (si) = 1, 
T s G 3 (s 1 ) = 0, 

|g 2 (s 2 ) - 0. 
|g 3 ( S2 ) = 0 . 

Gj(si) = 0, 

G<=(S2) = 1, 

GjfCsi) = 0 , 

G C 4 (S 2 ) = 0, 

|g 4 (si) = 0 , 



( 11 ) 


and that they span Pf[a] and P^[s] respectively. Accordingly we set the composite Lagrange-Hermite shape 
functions, for (a, s) e C, as (see also [3]): 


H^(a, s) = L c x (a) Gf(s), 
H^(a, s) = L c { (a) G^(s), 
H^(a, s) = L c { (a) G^(s), 
H%(a, s) = L c x {a)G c A (s), 


H^(a,s) - L%(a)G c x (s), 
s) = L^(a)G^(s), 
H^(a, s) = L^(a)G^(s), 
H^(a,s) = Lf(a) G^(s), 


and for (a, s) i C as H^(a, s) = 0, i - 1 ... 8. Refer to Fig. 3 for a graphical representation presentation of (12). 

It is easy to verify the following Lemma: 

Lemma 3.1. Let C given by (7) and 'Vc by (9). The following statements hold: 

• Every element v e 'Vc is uniquely determined by the eight Degrees Of Freedom (DOFs): point values v(V,) and 

s-derivatives f^v(Vj) at the four vertices V,-, i — 1 ... 4 of C. 

• The Lagrange-Hermite shape functions (12) constitute a basis for 'Vc- 

We are able now to define the interpolation over the cell C. In particular, let v e 'Vc and c,, d, the DOFs corre¬ 
sponding to the point v(V,) and ^--derivative values jfv(Vj) of v over the vertices of C. The scalar local interpola¬ 
tion function I c (-, •) 6 'Vc reads as: 


4 

I c (a, s) = ^ Cj Hfj_fa, s) + dH^fa, s), (a,s) e C . (13) 

i=i 

Similarly, we define the scalar global interpolation function, as a concatenation of the local interpolation functions 
I c (-, •), for all C e Ta^as- To that end, we first note that every internal node (a,, Sj ) of the discretization is a vertex 
to four computational cells, namely: 

Cj-lj-i = Oi-l,O',) X [Sj-uSj), Q-lj = Oi-i.a,) x [Sj, s j+ 1 ), | 

Cij-i - [ai,a i+ i) x [sj-i,Sj), Cjj = [a h a i+ i) x [s jt s j+l ). j 

It is moreover equipped with two degrees of freedom, v,j and d/j that correspond to the point and the .s--derivative 
values of the interpolated function respectively. In view of Lemma 3.1 the interpolation function can be uniquely 
(re-)constructed by the “full” set of DOFs. 
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Definition 3.1 (Scalar global interpolation function). Let 


{v,j, dij, i = 1... N a , j = 1... N s ] 

be a set of DOFs corresponding to point values and ^--derivatives over all discretization nodes. The interpolation 
function to these values reads as: 


/(«', s) = £ £ Vij (H^- 1 + H C 5 -'' j + H? J ~' + H C ; iJ ) (a, s) 

i=l 7=1 

+ ZZ ^ + H 6~ lJ + H A iJ " + H 2 J ) («. *) - 


(15) 


•=1 7=1 


for (a, s) e B. 

Similarly, the functions that belong piecewise in *Vc constitute the components of the scalar global approximation 
space: 


<V = 


B 


R 


v lc e 


VC ef, 


Aq'.A^ 


(16) 


We proceed to the 2 dimensional case and first introduce the 2D Lagrange-Hermite polynomial space as: 

V* = j(i ? x , v y ) T with v. v , v v e *y c | , 


(17) 


for C given by (7) and 'Vc by (9). We note that is a linear space of dimension 16. 

Given now the sixteen DOFs cf y , df y , i = 1... 4 that correspond to the point and s-derivative values of a function 
in 'V^ 1 the 2D local interpolation function reads as follows: 


4 


I c (a, s) = 


Z ( c ' H Z-i + d ‘ (a ’ 5) > 2 «^ i + < («. 

i=l i= 1 

4 4 

1 («>«) + ■ 


(18) 


1=1 


i=l 


Accordingly, follows the 2D global interpolation function, 

Definition 3.2 (2D global interpolation function). We define the vector local interpolation function as 

I(a, s) = (f x (a, s), I y (a, s)) T (19) 

where I xy are the corresponding scalar global interpolation functions. 

After invoking (15) and the global point and ^-derivative values d X ’J, i = 1... N a , j = 1... A v }, (19) reads 

as 


N„ N s 


i= 1 7=1 


'■=1 7=1 


+ Hy J -' + H C f ’) (a, s) 
-//' -//; )«r. v). 


( 20 ) 


The corresponding 2D global approximation space reads: 

= jv : B —4 R 2 1 v| c e VC e r A „, Aj } . 


( 21 ) 


Based on the strong formulation (3) and on (21) we present here, without further justification, the corresponding 
FE formulation for the numerical treatment of (3) and refer to [14] for further details on the “full” FEM. The 
presentation is with respect to the family F; symmetrically similar is the problem for the other filament family. 
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Find continuous in time solutions F(-, •, t) e < V 2d such that, for every test function v 6 'V 2d , 
0 = f d(p B d 2 s F ' + p A D,F ■ v + /WAF ■ <3.,v) d(a, s) 

B 

+ f - 0o)3,s-F x • 8 s y + p s (D ,F - D,*F*) ■ v) d(a, s) 

- J p(Q){d a ^ ■ d s \ - ■ <9„v) d(a, s) 

+ f d(f an(a)d s F + /inn(a)v ± (a)) • v da 

JdBn\s=-L} 

*r 

JdB 


JdBn{s=-L) 

^tether ' V da . 


( 22 ) 




3.2 Numerical considerations 


In this section we present some of the considerations that we took into account when addressing the numerical 
resolution of (22). In particular we address the numerical treatment of each of the terms of (22) separately. 


Following (18), F is represented in terms of the Lagrange-Hermite basis functions (12) over the computational 
cell C 6 Taq- as at the time step t k as 


F k \ (a,s) 
I c 



(23) 


where p k ix , p k y i — 1... 8 denote the x and y respectively, point and derivative value DOFs assigned to the vertices 
of C under lexicographic enumeration, see also (8). There are some notable exceptions in this approach where it 
was deemed necessary to approximate the solution with the help of another basis set. We will address these terms 
with detail. As a rule though, for the derivation of the FEM, we replace (23) in (22) and test against every shape 
function Hf for k - 1... 8 and C 6 T\ llt ,\ s . 


3.2.1 Filament length 


We compute the length of the filaments L, and the length distribution tj explicitly, following the modelling consid¬ 
erations described in [15] (see also [5]). The filament length L, in particular, is given by 


L 


^cap,eff 

^sev 


^cap,eff 2v 
-T— + - 


n(s = 0) 


where K sev , /c C a P ,eff are the severing and capping rates of the filaments, and where ? 7 m ; n is a minimum cut-off density 
for F-actin. Note moreover that faster polymerization rate v leads to wider lamellipodia. 


3.2.2 Filament bending 


The bending term reads in the FE formulation (22) as 



■ d 2 s \ d(a, s) . 


The numerical solution F is disretized implicitly in time and is represented over the computational cell C at 
the time step r' !+1 using (23). It’s contribution in the FEM formulation reads —after testing against H^(a,s), 
k = 1... 8— in its x-coordinate as follows: 

^IX +1 


Li 


diH ( r d]H^ 


) (a, s) d(a, s) . 
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3.2.3 Adhesion with the substrate 


The adhesion term reads in the FE formulation (22) as 


X 


tjD ,F • v d(a, s) , 


which reads after explicit in time discretization, as: 

( JT W+ 1 — JT n 


X 


At 


vd s F" ■ v d(a , s ) . 


Expanding F", F" +1 according to (23) and testing against //', it reads in the ^-coordinate as 

8 8 

£(£j^ (^' _ p '-<) H > (a ’ s) ~ v X Pi* 9 * 11 ?' (a ' ^) H k (a ’ s) d(a ’ s) ’ 


or 

1 

At 


8 8 

J pfc 1 £ ( H f H k) (<*. *) S) - J Pl(£ t £ (flf ) (a, *) rf(ar, *) + V £ (d s H? H c k ) (a, s ) d(ar, s)). 


The implicit terms are added in the stiffness matrix, and the explicit ones in the right-hand side of the discrete 
FEM formulation. 


3.2.4 Stretching of cross-links 

The stretching of cross-links reads in the FE formulation (22) as 

£ r/rf (D,F - D*F*) ■ v d(a, s) . (24) 

It involves filaments of both families and is derived under the assumption, that in a previous time t—6t the filaments 
a, a* cross at the positions s, s*. 

Due to its nature, (24) necessitates special treatment: for every filament and discretization node (a, s) of one 
family, we assume that there exist (a*, s*) of the other family (in general non-discretization nodes) such that: 

F"(a, s + vA t) = F"’*(a*, s* + v*A t) , (25) 

where F" denotes the numerical solution of the F at time t". Accordingly, (24) reads, after time discretization, as: 

Af(AF(a, s)-D,F*(tt*,j*)) 

=F" +1 (of, s) - F"(or, s) - vAtd s F n (a, s) 

- F" +I, *(a*, s*) + F"'*(a*, /) + v* Atd s F n '\a\ /) 

=F" +1 (a, s)-F"(a,s + vAt) 

- F" +1 ’*(a*, s') + F”’*(or*, 5* + v*At) 

=F" +1 (or, s) — F" +1 ’*(or*, s*) . (26) 

The relation (26), necessitates that for every (a, ,v)-pair (discretization nodes) the corresponding (a*, ,v*)-pair (dec¬ 
imal non-discretization nodes in general) should be computed. We do this in the following way: we approximate 
the crossing-point in terms of the (*) family, by identifying the “below” and “above” (bounding) filaments and 
position nodes a* { , a* 2 and .y*, s* respectively, see also Fig. 2. These provide with the four surrounding vertices of 
the crossing-point with respect to the (*) family: 


(<7i.o<?iv) T = F*(a*,s*), 
(tttx, <73>’) t = F*(a*,s*), 
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(d2x, dlyV = F*(a*,S*), 
(q4x, d4yV = F*(a*,S*). 



The corresponding point is approximated by the interpolation: 

F‘(a*,4)*(^,^) T 

= CiF*(aq, jj) + c 2 F*(ai, 4) + C3F*(«2> 4) + c 4 F*(ar 2 , * 2 ) , 

with the weights given through the integer parts (a, s ) of (a*, s*)\ 

Ci = (1 - a)( 1 - s), c 2 = (1 - a)s, 

C 3 = or(l - s), C 4 = as. 


(27) 

(28) 


(29) 


Subsequently, we multiply with the Lagrange-Hermite shape functions Il,(a. s), i = 1,..., 8 (12) defined over the 
cells of the domain B of F. Since though the nodes of one family are represented only as decimal coordinates with 
respect to the other family, one rectangular cell in B corresponds to a quadrilateral in B* and covers there parts 
of different cells. This causes the following complication: if the “usual” Lagrange-Hermite elements are used 
for the (*) family in the reconstruction of F* over B, it becomes unclear which direction should the derivatives 
take into account. It also turns out that the use of the more “usual” Lagrange-Hermite approximation for the F is 
problematic because the two parts are then implemented in an asymmetric way (the term acts of the points and the 
derivatives of F, but only on the points of the (*) family). We have used instead Lagrange-Lagrange 1 elements for 
the approximation of both families. 


Therefore, instead of using (23), we expand F and F* in the cell C of B as: 

4 4 3 t 


F|^(a, s) = 


F*| c (a, s) = 


1=1 1=1 

4 _ 4 _ 

J]q»; x l Lf(a,s), ^^ + 1 Lf(a,4 


V 1=1 


i= 1 


n n +1 «.n+l „n +1 „«+1 


q“y are the DOFs corresponding to the positions of the vertices of C with respect to the 

—n+l _ v4 n +1 „n +1 r .^A — M +1 


where p^, q[ x , _ ±iy 

two families. Due to the interpolation (28) we write q"* 1 = 2r=i tfVtflrx anc * ' = E?=i w > t h q"^, 

q" + r ' , c n A 1 , r = 1... 4 given by (27) and (29) respectively. 


At the end, after testing against Hk, and integrating over C, (26) reads in the ^-coordinate as: 


J (F" +1 


• F " +1 ’*)H k 


z 


■ y c n+i a n+i 

/ i ^ i,r Hi,rx 


J' (if Hk) (a, s) d(a, s) . 


3.2.5 Twisting of cross-links 

The twisting of the cross-links reads in the FE formulation (22) as 

J" qq\4> - 00)3.^ ■ d s \ d(a, s) . 

The angle <p between the crossing filaments is approximated with a process similar to the one described in Sect. 
3.2.4. At every discretization node (a, s) of one family we identify the (probably decimal) node (a*, s*) of the 
other family such that F(cr, s) = F*(a*, s*). We employ the interpolation formulas (28) adjusted for the d s : 

d s F*(a\ s*) =cid s F*(a*, s*) + c 2 3 s F>*, s* 2 ) 

+ c 3 d s F*(a* 2 , 4) + c 4 3 J F*(ar*, 4), (30) 

with weights ci... c 4 given by (29). Expanding now f| implicitly (at t" +1 ) according to (23), and computing cf> 
explicitly (at t "), the twisting term reads in its ^-coordinate as 

8 

-m* ~ Mp'Cy 1 

i=l 

1 linear shape functions in both a- and ^-directions of the form L^(a) L c .{s), i,j= 1,2 with L c (-) given by (10) 


J r (3,v//f d s H ^) (a, s) d(a, s) . 
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3.2.6 Filament repulsion 


The pressure term reads in the FE formulation (22) as 

J' p(e) (dffF - 1 • d s \ - d s F ± ■ d a \j d(a, s ) 


where 


Q = 


13,F • d a F±\ 

pie) = p p e ■ 


(31) 

(32) 


For the computation of g , and in effect of pig), we approximate the point values F(cr, s)| in (31), by the —explicit 
in time— cell averages 


1 


AsAcr J c 


g 

£ F "(a, s) d(a, s) = ^(p" r ,p" v ) T £ Hf(a, s) d(a, s ) . 


(33) 


After implicit discretization and expansion of F according to (23), the pressure term reads in the x-coordinate 


as: 


- 2 pVv'p^ (£ { d « H f d ' H k ) («. s) s) + £ (3^, C d a H f) (a, s) d(a, s)J , 


with p(g) given by (32), (31), (33). 


3.2.7 Innerforces 


The computation of the inner pulling forces f m (a) and follow from the force balance law (see also the BCs 

(5) at s = -L): 


f r](a, s)(f tan (a)d s F(a, s) + f inn (a)Y(a)) 
J\0,2 tt) 


)[0,2 it) 

In effect, this leads in a variational approach to the minimization of 


da - 0 . 


s——L 


f rj(a , s)(- (/tan (ff) - 7A) 2 + (/inn (a) - (1 - y)A) 2 ) 

J[0,2n) 1-7 ’ 


da 


s——L 


(34) 


under the constraint (34), where A = g 1F (A c - Ao)+ measures the positive deviation of the area A c occupied by the 
cell from an equilibrium value Aq. The constraint minimization problem yields 


/an(o') = 7A 1 1 - ipi,p 2 ) T ■ d s F(a, s) 
/inn (or) = 1 - yA(l - (pi,p 2 ) T ■ V(q')) , 


s——L 


where 


(pup 2 ) T = 


f . 

0,2tt) 

</ ■ 


q (y<9. s F ® d s F + (1 - y)V(ar) ® V(a)) 


da 


s=-L 


-1 


qiydsF + (1 - y)V(ar)) 


da 


where (x) stands for the regular number multiplication. 

Using (35) and sq = -L, we treat the first term of (34) as follows (the second term is treated similarly) 
I tj (a, s 0 )yA [d s F(a, so) - (p\,P 2 ) T ■ d s F(a, s 0 )3,F(q;, s 0 )j da , 

J[ 0 , 2 k ) 


(35) 

(36) 


( 37 ) 
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or, after expanding F 


I c 


explicitly at t" according to (23), testing against //;., it reads in the x-coordinate as 


qyA 


Y p l* f ( d s H i H k)(a,so ) da - Y pIx({pIx’ p %) ' 0“i,A*2 ) T J J'(d s Hi d s HjH k )(a, s 0 ) da 


\i= 1,..,8 


where the integrations are over each inner element of the discretization and where (fj\ ,/T 2 ) t is computed explicitly 
in time via (37). 


3.2.8 In-extensibility 


The in-extensibility term reads in the FEM (22) as 


J' r)A inexAF • d s \ d(a, s) . 


We employ the Augmented Lagrangian approach to evaluate /lj next ; the strong formulation of which recasts into 



with A" been updated in every time step by 

A” +1 = A" + - (<9 S F" • <9 S F" +1 - l) , 


for 0 < e < 1. When linearizing as <9jF i—-> (|3.,F| 2 - l) <9,F at 3,F'', the punishing term -^<9^ (|<9jF| 2 - l) <9 s Fj is 
approximated by 

-^((|3 S F"| 2 - 1) <3, S F" +I +2(3, S .F" • ^F' 1+1 - |3 S F"| 2 )5 S F"), 
revealing a contribution in the (indirection of both the old and the new time step. 

In effect, the in-extensibility term reads as: 

-d s J(,l" + (|3,F' ! | 2 - l) )(9 S F" +1 + ? (3jF" ■ 3 S F' !+1 - |<9 S F"| 2 ) 3 S F"J , (38) 

with A" given by 

A n+] = A" + -(|3 S F"| 2 - 1) + ~(d s F" ■ 8 s F n+l - |<9 S F"| 2 ) . (39) 


In a way similar to the previous terms, we expand F"| , F" +1 | according to (23), A" by £r=i 'l™. and we test 
against Hu to obtain the final form of the in-extensibility term (38). The formula contains both explicit and 
implicit contributions and is omitted for the sake of the presentation. 


3.3 Reparametrization 

The direct implementation of the FE formulation (22) as was described in the previous paragraphs has been seen 
to have two drawbacks that we have addressed with proper reparametrizations, at every time step, along the a- and 
i-directions separately. We sketch here the suggested reparametrization and the treatment that we have provided 
and refer to [14] for further detail. 
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3.3.1 In a-direction 


We reparametrize along the a-direction to ensure that the “computational” filaments (mappings of the discretiza¬ 
tions nodes along the a-direction) are “regularly” distributed over the physical space il(t) (see Fig. 1). In particu¬ 
lar, we define a mapping g(-, t) : [0, In) —> [0, 2n), with (5 \-> a := g(/3, t), which leads to a new weak formulation 
of the problem. The original bending and adhesion terms, for example: 

IUfsd 2 s (r]d 2 s F) + fi A T]D,F 


recast into: 


where F(a, s ) = F(/3, s), and 


IUfid 2 (f]d 2 F) + i^tjDPf 

s ) = T](a, s) |g'(£)| , 
D) - d, - vd s , 


(40) 

(41) 


where (41) is derived under the simplification assumption that g 1 « 0. Note also that, in (40), we have incorpo¬ 
rated the contribution of g into f/. 


One way to define such g that maintains the “regular” distribution of the “computational” filaments is the follow¬ 
ing: 

[J \d a F(&, 0, t)\da ' 


a i * g (a, t) =: p ■ 


\d a F(a, 0, t)\ da 


2 n , 


which for M(t) = f Q \d a F(a, 0, f)| da reads | d p F(p, 0,01 = ^. 


Numerically, this is achieved by setting Li = 2n, and defining the piecewise linear g: 


8 ■ 


j= 1 7=1 


[(/ - 1) Aa, i Aa] 


P i—» a := Aa(/- 1) + Aa 

where Aa is the discretization step size along the a-direction. 




3.3.2 In v-direction 


Uneven filaments lengths, lead to discretizations with uneven As in the opposite sides ai, ai of the cell C, see 
(7). In some cases this can lead to numerical instabilities in the form of filament oscillations. To address this 
issue we have remapped all filaments to a constant length and transferred the control of the filament length to the 
in-extensibility terms. 


In particular, we set a map from [-L(a, t), 0] onto [-L, 0] for a constant L via the change of variables: 

(s, a) - 


L(a) 


, a 


where L(a) is the ratio between the original constant length L and the new length. Accordingly, the strong formu¬ 
lation of the problem (3) is transformed to read 


d 2 s (r,d 2 s F) + (L(a))*Jd t --^- ) d)F 

+ (L(a)) 2 d s (r ir] *(c!>-<Po)d s F^) 

+ " "• <««» 4 (( ft ~ik d -) r -(*'-zM r ) 

- d s (r]A inext d s F) 

+ (L(a)) 3 (5, (p<9 0 ,F J_ ) - 8 a (pi5jF J_ )) = 0 . (42) 
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Figure 4: Movement of a cell on an adhesive substrate (red) with less-adhesive stripes (white). Shading represents actin 
network density. A-l: Times series of a cell moving over a stripe pattern (80% drop in adhesiveness). Parameter values as 
in Table 1. The bar represents 10 pm. 


where the (punishing) in-extensibility ,lj next now reads 

dinext = ^ (|3 f F| 2 - ( L(a)f ), 

accounting hence for the proper length of the filaments. This term is treated numerically in the way described in 
Sect. 3.2.8. The appearing in the stretching and adhesion terms is absorbed by redefining the polymerization 
velocities. Note also that L(a) does not depend on s, hence all integrations by parts in the s-direction are not 
affected. 


4 Numerical simulations 

The purpose of this section is to demonstrate that the model is capable of describing complex biological experi¬ 
ments. In [15] several numerical experiments were described which demonstrated the effects of the the individual 
terms of the model. Additionally we demonstrated that the model can be used to simulate chemotaxis and to study 
how different signaling processes affect cell shape and filament density. Here we go a step further and simulate 
how the shape of a migrating cell is influenced by spatially selective adhesion patterns. Such studies are used to 
better understand the force balance between adhesions, contraction mechanisms, actin polymerization and other 
network proteins. 


4.1 Experiment 1: Adhesive\Less-adhesive stripes 

In this series of experiments we show that the model’s behavior is similar to that of cells in the experiments de¬ 
scribed in [4]. In these experiments migrating fish keratocytes were placed on substrates which were prepared to 
have chemical patterns of the Extracellular Matrix (ECM) protein, fibronectin. This protein is a ligand of integrin, 
an important component of adhesions. In [4] stripe pattern were used with 5 pm wide adhesive (fibronectin con¬ 
taining) stripes and non-adhesive stripes (without fibronectin) with widths varying between 5 pm and 30 pm. In 
[4] it was described, that cell shape was affected in a very distinct way: protruding bumps on the adhesive stripes 
and lagging bumps on the non-adhesive stripes were observed and their width was correlated to the stripe width. 
Also it was observed that cells tend to symmetrize themselves such that they had an equal number of adhesive 
stripes to the right and to the left of their cell center. In the numerical experiment we used the same geometrical 
pattern with 5 pm adhesive stripes interspaced with 7 pm less-adhesive stripes. In the mathematical model, the ad¬ 
hesion forces result in friction with the ground and, speaking in numerical terms, they link one time step with the 
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Figure 5: Movement of a cell on an adhesive substrate (red) with less-adhesive stripes (white). Shading represents actin 
network density. A: 90% drop in adhesiveness, B: 80% drop in adhesiveness. Parameter values as in Table 1. The bar 
represents 10 |im. 


next. Therefore we cannot locally set the adhesion coefficient to zero, since this would render the stiffness matrix 
degenerate. Hence the adhesion coefficient in the less-adhesive regions was decreased to 10 - 20% as compared 
to the adhesive regions. Whilst the keratocytes in the experiments by [4] move spontaneously without an external 
signal, we simulate chemotactic cells, since at this point the model cannot describe the dynamics of contractile 
bundles observed in keratocytes. However the numerical results show, that there are many similarities as far as 
general behavior and morphology are concerned, suggesting that the underlying phenomena are very similar. In 
Fig. 4A-I a time series of the behavior of a cell on stripes with a drop in adhesiveness of 80% is depicted. The 
following agreements between the simulation and the experiments were found: 

• On the striped region the cell shaped became more rectangular as compared the the crescent shape in the 
homogeneously adhesive region. 

• The numerical cell shows protruding bumps on the adhesive stripes and lagging bumps on non-adhesive 
stripes. 

• The simulation started with the cell being slightly non-symmetric with respect to the adhesive stripes and 
changed its position to have an equal number of adhesive stripes to the right and to the left of its center. 

• Spikes at the cell rear were observed. 

• After leaving the striped region the cell resumed its crescent shape and continued to migrate as before. 

In Fig. 5A and B a comparison between the bumps on stripes with a 90% (A) and a 80% (B) drop in adhesiveness is 
shown. Here the ^--discretization used was twice as fine to allow for better analysis. As expected the bumps of the 
90%-drop stripes are more pronounced. Over a time interval of several minutes we also observed the fluctuations 
in bump width observed in [4], 


4.2 Experiment 2: Less-adhesive spikes on strongly adhesive ground 

In the next simulation we describe the behavior of a migrating cell on a substrate with a different pattern, showing 
that the model can make predictions for future biological experiments. We use the same setup as above with a 
pattern which consists of two shifted less-adhesive spikes from above and below. The drop in adhesiveness was 
chosen to be 80%. As opposed to the situation above, the cell is now able to almost fully avoid the less-adhesive 
regions. The behavior observed over a time of 30 min is depicted in the time series in Fig. 6. It can be seen that the 
cell behaves as if the less-adhesive spikes were obstacles and always only a very small fraction the lamellipodial 
region enters the less-adhesive areas. 
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Figure 6: Movement of a cell on an adhesive substrate (red) with less-adhesive spikes (white). Shading represents actin 
network density. Parameter values as in Table 1. The bar represents 10 pm. 


4.3 Parameters values: 

For the discretization we used a time step of 0.002 and nine nodes per filament. For the first experiment we used 
36 and 72 discrete filaments, for the second one 36. 

For the biological parameters, we used the same as those in [15], apart from the adhesion coefficient which was 
increased for the adhesive regions and decreased for the less-adhesive regions. They are summarized in Table 1. 

Aknowledgements: N.Sfakianakis wishes to thank the Alexander von Humboldt Foundation and the Center of 
Computational Sciences (CSM) of Mainz for their support, and M. Lukacova for the fruitful discussions during 
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